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ABSTRACT 

The reionization history of the intergalactic medium (IGM) at high redshift (z Jj 6) was Ukely strongly shaped 
by several global feedback processes. Because the earliest ionizing sources formed at the locations of the 
rare density peaks, their spatial distribution was strongly clustered. Here we demonstrate that this clustering 
significantly boosts the impact of feedback processes operating at high redshift. We build a semi-analytical 
model to include feedback and clustering simultaneously, and apply this model to the suppression of star- 
formation in minihalos due to photoionization. The model is built on the excursion-set-based formalism of 
Furlanetto, Zaldarriaga, & Hernquist (2004), which incorporates the clustering of ionizing sources, and which 
we here extend to include suppression of star formation in minihalos. We find that clustering increases the 
mean HII bubble size by a factor of several, and it dramatically increases the fraction of minihalos that are 
suppressed, by a factor of up to ~ 60 relative to a randomly distributed population. This enhanced suppression 
can significantly reduce the electron scattering optical depth t, as required by the three-year data from the 
Wilkinson Microwave Anisotropy Probe (WMAP). We argue that source clustering is likely to similarly boost 
the importance of a variety of other feedback mechanisms. 
Subject headings: cosmology: theory - early Universe - galaxies: high-redshift - evolution 



1. INTRODUCTION 

The three-year data on the polarization of the cos- 
mic microwave background from the Wilkinson Microwave 
Anisotropy Probe (WMAP) satellite have provided an impor- 
tant new constraint on the reionization history of the inter- 
galactic medium (IGM). The WMAP experiment has mea- 
sured the optical depth r to free-electron scattering between 
us and the recombination epoch at z « 1100. The result, 
r = 0.09 ± 0.03, implies a redshift for instantaneous reion- 
ization of z « 11 (Spergel et al. 2006), and can be used to 
rule out a tail of significant ionization extending out to higher 
redshifts. Another constraint on reionization follows from the 
Gunn-Peterson (GP) troughs detected in the spectra of bright 
quasars at z « 6. These GP troughs probe the neutral hydro- 
gen along the line of sight, indicating that the last ^0.1% 
of the neutral fraction was rapidly disappearing around this 
epoch (Fan et al. 2006; Mesinger & Haiman 2004). 

The relatively short reionization history implied by the 
combination of these constraints requires a strong suppres- 
sion of ionizing radiation from "minihalos" with virial tem- 
peratures Tvi, < 10** K (Haiman & Bryan 2006). A variety of 
feedback mechanisms might play a role in this suppression 
(e.g. Haiman & Holder 2003). In this paper, we will explore 
one such feedback mechanism - the suppression of gas infall 
and star formation in minihalos in photoionized regions - and 
focus on how the clustering of the ionizing sources enhances 
this effect. 

Many other feedback mechanisms may also play a role in 
driving the history of reionization. Some (such as the buildup 
of an X-ray background) may affect the whole IGM more or 
less uniformly and simultaneously, but most will be localized 
to within a small fraction of the Hubble distance, at least for 
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some time. Examples of such localized feedback mechanisms 
are (1) the first generations of stars will enrich the intergalac- 
tic medium with metals (Madau et al. 2001; Aguirre et al. 
2001), which will affect the formation and evolution of subse- 
quent generations of stars (Tumlinson & Shull 2000; Bromm 
et al. 2001; Schaerer 2002; Wyithe & Loeb 2003); (2) ioniz- 
ing UV photons can both catalyze the formation of (Ricotti 
et al. 2002b; Oh & Haiman 2002; O'Shea et al. 2005), and di- 
rectly destroy (Haiman et al. 1997a) H2 molecules, which will 
affect the ability of gas to cool and form structures (Haiman 
et al. 1997b; Haiman et al. 2000; Ciardi et al. 2000), (3) gas in 
regions that have been ionized may not be able to settle into 
the shallow potential wells of small dark matter halos, exac- 
erbating the H2-photodissociating effects of a Lyman Werner 
background (Haiman et al. 1997b; Mesinger et al. 2006) and 
increasing the minimum size required for a halo to contain 
collapsed baryonic structures (Ricotti et al. 2002a,b; Oh & 
Haiman 2003). 

Since these mechanisms operate over a limited length scale, 
their effects will depend strongly on the spatial distribution of 
halos hosting ionizing sources. Numerical simulations are a 
promising way to address feedback among clustered sources, 
since they capture the full, three-dimensional relationships 
among the host halos. However, the dynamic range required 
to resolve the small minihalos, within a large enough cosmic 
volume to be representative, remains a challenge, especially 
in simulations that include radiative transfer (see, e.g. Iliev 
et al. 2006). Semi-analytical models avoid the problems as- 
sociated with the large dynamical range; they are also an effi- 
cient way to explore parameter space and serve as important 
sanity checks for more complicated simulations. 

We believe that this paper describes the first semi-analytical 
model to include feedback and clustering simultaneously as a 
way to suppress ionizing radiation from minihalos. Prior work 
on the clustering of halos has focused on the problem of de- 
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riving analytical expressions for the bias between the distribu- 
tion of dark matter halos and of mass (Scannapieco & Barkana 
2002); (see also Mo & White 1996; Scannapieco & Thacker 
2005). These results have been applied to the distribution of 
minihalos at high redshifts: Haiman et al. (2001) included the 
bias to account for the effects of minihalo photoevaporation; 
Iliev et al. (2005) likewise characterized the effect of minihalo 
clustering on the propagation of ionization fronts through the 
IGM. Minihalos were considered to be sinks of ionizing pho- 
tons in both of these works (as well as in Ciardi et al. 2006; 
Furlanetto & Oh 2005). In contrast, in the present paper, we 
examine minihalos as potential sources of ionizing radiation. 
This means that the radius of influence of our feedback mech- 
anism cannot be simply evaluated at a fixed time-step; instead, 
it is a function of the entire ionization history of the region, 
including the history of mergers between ionized regions. The 
preceding methods cannot be easily be adapted to such a cal- 
culation, and we have been driven to employ a novel Monte- 
Carlo means of simulating the coupled history of ionization 
and feedback around any given point. Our primary conclu- 
sion is that, because of their highly clustered distribution, lo- 
cal feedback effects can dramatically suppress star formation 
in minihalos. 

The rest of this paper is organized as follows. In § 2, we 
describe our model for reionization, extending the formalism 
of Furlanetto, Zaldarriaga, & Hernquist (2004) to include an 
approximate treatment of biased halo distribution and feed- 
back. In § 3, we report our results for the reionization history, 
and discuss the impact of clustering on the overall suppres- 
sion of minihalos, on the size distribution of ionized bubbles, 
and on the optical depth to electron scattering. In § 4, we 
discuss several limitations and possible extensions and future 
applications of our model. Finally, in § 5, we summarize our 
conclusions and the implications of this work. 

2. MODELING REIONIZATION WITH CLUSTERING AND 
FEEDBACK 

We have chosen to implement the suppression of minihalos 
in ionized regions (Thoul & Weinberg 1996; Dijkstra et al. 
2004; Oh & Haiman 2003), both to understand the effects 
of this particular mechanism, and as an example of localized 
feedback mechanisms in general. We parameterize the sup- 
pression effect by setting different minimum virial tempera- 
tures for a halo to contribute to reionization rvir.min, in ionized 
and neutral regions: 



10 K in neutral regions, 
10^ K in ionized regions. 



(1) 



We will hereafter refer to halos with virial temperatures be- 
tween 10^ K and 10'* K as "minihalos", and halos with virial 
temperatures above 10"* K as "large halos". 

Furlanetto, Zaldarriaga, & Hernquist (2004) (hereafter 
FZH) have created an elegant, excursion-set-based approach 
to following the ionized fraction of the Universe, which cap- 
tures the clustering of ionizing sources. They begin with the 
assumption that the amount of mass, Mion, that sources within 
a given halo can ionize, is simply proportional to the total 
mass of the halo, 

Mion = CMhdo, (2) 

where ( is the "ionization efficiency factor". The ionized frac- 
tion x; of a region is then proportional to the fraction of mass 
in that region that is contained in collapsed objects, 

Xi = CfcowiM, (5,Z, Tvir.min), (3) 



where M is the mass of the region, S = (p- p)/p is the over- 
density of the region, and z is the redshift. 

The collapsed fraction in the extended Press-Schechter for- 
malism is: 



/coii(M,(5,z,Mmi„) = erfc 



S,(z)-S(M) 



v/2[a2(M^i„)-a2(M)] 



(4) 



The critical density for collapse is (5e = 1 .686/g(z), where g(z) 
is the linear growth factor of density perturbations, normal- 
ized to unity at z = (we use the approximation from Carroll 
et al. 1992). The variance of d(M), cP-{M), is determined from 
the power spectrum, which we in turn calculate using the fit- 
ting formulae by Eisenstein & Hu (1999). 

We adopt the relation between virial temperature and halo 
mass given by Barkana & Loeb (2001): 



r^i^(M,z)=1.98xlO^K 



M 



n. 



A, 



lOi^/z-iMo 

1/3 



2/3 



1+Z 

10 



riMiz) 187r2 



0.6 



(5) 



where Ac EE 1 Stt^ + 82[ 1 - f7M(z)] - 39[ 1 - f7M(z)]^ and we set 
^ = 1.22 for Tvir < 10^ K and M = 0.59 for T^i, > 10"* K. We 
solve this expression for M to determine the minimum halo 
mass, Mjnm(Tvh.min,z), as a function of redshift. 

The condition for a region of mass M to be fully ionized 
is C'/coU > 1- FZH solve this equation for S{M) yield- 
ing a barrier function SxiM,z)- The largest mass at which 
S(M) > Sx(M,z) about a particular point is then identified as 
the mass of the ionized region to which the point belongs. 
Computing the mass M defined by this identification at dif- 
ferent redshifts then translates into the history for the mass 
of the ionized bubble in which the given point resides. Con- 
versely, the mass function of ionized bubbles at any redshift 
can be constructed either by Monte Carlo realizations of many 
S(M) trajectories, or from the statistics of the first-crossing 
distribution of a random walk (for a method of finding the 
first-crossing distribution for an arbitrarily shaped barrier, see 
Zhang & Hui 2005). 

The cosmological parameters used throughout this paper 
are: h = 0.72, flufi = 0.27, f^A = 1 - ^Mfl, f^b.o = 0.045, 
(78 = 0.9, n = 1, taken from the best-fit values in the 1-year 
WMAP data release (Spergel et al. 2003). For our purposes, 
the most significant changes in the cosmological parameters 
in the 3-year WMAP results (Spergel et al. 2006) are the lower 
values of crs = 0.74!|j ;J]^ and n, = 0.95 l^g;;]!!. The main effect 
of both of these changes is to slightly delay structure forma- 
tion and reionization. This can be compensated for by in- 
creasing the efficiency parameter C,, and it does not affect our 
primary conclusions about the effect of clustering. 

2.1. Adding feedback 

We next depart from the FZH formalism to incorporate pho- 
toionization feedback. Adding feedback means that the bar- 
rier S(M) at any given redshift can no longer be pre-specified; 
rather, it will now be a function of the entire ionization his- 
tory of the region being considered. We therefore must resort 
to Monte-Carlo simulations to calculate the coupled history 
of ionization and feedback around any given point. 

If we assume that minihalos are randomly distributed and 
unbiased with respect to the ionized bubbles, then the ionized 
fraction of a region can be found by integrating the following 
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Fig. 1 . — The ionization history of the Universe: mean ionized fraction versus redshift. The solid curve is the history with biased feedback. The leftmost 
short-dashed curve (blue in the color version) includes contributions from all halos, the rightmost short-dashed curve (red) includes only large halos, and the 
upper long-dashed curve (purple) includes unbiased feedback. All of the above use C, = 12. The light-weight long-dashed curve (purple) is an unbiased feedback 
calculation with a lowered efficiency factor, C = 5. Comparing the lower unbiased feedback curve with the biased feedback curve demonstrates that minihalo 
clustering must have a large effect on the early stages of reionization even if the reionization redshift is assumed to be known. 



simple differential equation, starting from high redshift: 

dxi 
"5F 



"■^i, large ,, > "-^i.mini 

■+(1-Xi)- 



dz 



dz 



(6) 



where c/xiia,ge/iiz is the contribution to ionizing the region 
from large halos and dxi„iim /dz is the contribution from mini- 
halos: 



'''''''"^'=-C%i(rv,->io4K) 



dz 

"■^i,mi 



dz 



(7) 



dz 



^^:^(102k < Tvir < lO'* K) (8) 

dz 

Note that these expressions follow FZH and ignore recombi- 
nations (we will discuss this issue in § 4.1 below). Applying 
these equations to the whole Universe (i.e., to a region of mass 
M ^ oo, 6 = 0) gives us the reionization history of the Uni- 
verse in the absence of any clustering. 

The efficiency factor is chosen to be C = 12, unless stated 
otherwise. This factor can be thought of as the product 



C = /esc/*A^7/i (1 + "rec) 



(9) 



where /esc is the fraction of ionizing radiation that escapes 
from halos, /, is the fraction of baryons that end up in stars, 
Ny/i, is the number of ionizing photons produced per stellar 
baryon, and Wrec is the average number of times a hydrogen 
atom recombines. As an example, (^ = 12 could be obtained by 
adopting the reasonable choices of /esc = 0.1, /* = 0.2, N^/i, = 
4200, Hiec = 6. We ignore the dependence of Wiec on bubble 
size, an approximation which becomes inaccurate toward the 
tail end of reionization (Furlanetto & Oh 2005). We discuss 
this issue further in §4. 1 . 

The integration is performed numerically from z = 35 to 
z = 7 with step size dz = 0.02. All distances and densities 
specified in this paper are in Lagrangian comoving coordi- 
nates. Lagrangian coordinates are appropriate for our prob- 
lem since we care more about the column density between 
two points (which is directly related the required number of 
ionizing photons) than about their physical separation. 

Figure 1 shows the evolution of the ionized fraction in the 
above toy model. The figure displays the history with only 
large halos allowed to contribute (rightmost short-dashed. 



red, curve), with all (i.e. both large and mini-) halos con- 
tributing (leftmost short-dashed, blue, curve), and with mini- 
halos suppressed by unbiased feedback (upper long-dashed, 
purple, curve). The first two curves will bound the reion- 
ization history with feedback, since they represent the most 
extreme cases. As expected, the unbiased feedback history 
begins following the all-halo curve, then lags behind as the 
ionized fraction increases and minihalos are suppressed. 

2.2. Capturing clustering 

In order to capture the biased distribution of the minihalos 
with respect to the ionized bubbles, we begin by synthesizing 
random walk trajectories representing the overdensity S{M) 
around a random point in space averaged over spheres of mass 
M. Figure 2a shows such profiles. The overdensity tends to- 
wards zero as M gets large, since density fluctuations average 
out in very large spheres. 

We then divide the trajectory up into thin spherical shells 
and numerically integrate equation (6) for each shell sepa- 
rately. In this paper shells are logarithmically spaced from 
lO'*-lOi**M0 with thickness £/log(M/Mo) = 0.015. Shells 
that overproduce ionizing photons are allowed to contribute to 
the ionization of the outside neighboring shell (this is calcu- 
lated recursively to allow extra photons to spill outward until 
consumed). For the purpose of finding the mass of the ion- 
ized bubble surrounding the central point, as defined above, 
the ionized mass is still associated with the shell from which 
the ionizing photons originated. 

The total mass ionized by sources within a sphere of size M 
is the sum of the contribution from all interior shells. 



M,o„(M) = J2 ^.on,shell(M'). 



(10) 



M'<M 



As in the FZH formalism, a point is said to belong to an ion- 
ized bubble of mass M such that M is the largest mass for 
which Mion(M)>M. 

Figure 2b shows the mass histories corresponding to the 
four density profiles shown in Figure 2a for the all-halo, 
large-halo-only, and biased feedback cases. As was found 
in Furlanetto & Oh (2005), the growth of the region proceeds 
primarily by large jumps in mass, corresponding to "major 
mergers", wherein a point is engulfed by a neighboring, larger 
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Fig. 2. — Examples of the density profile (a) and the corresponding HII region mass history (b) about four randomly chosen points. Heavy (black) lines in (a) 
show the mean overdensity in a spherical region versus the mass of the region, extrapolated linearly to the present. Light (red in the color version) curves mark 
the itlcr level for the density fluctuations. In (b) the HII region mass to which each point belongs is plotted against redshift for three cases: (left) all (mini- and 
large) halos contribute to reionization, (center) feedback is turned on, and (right) only large halos contribute. 



HII region. The mass of a merger product is determined by 
where a peak occurs in the &{M) trajectory, so many of these 
jumps are repeated in each history at the same mass, but at 
different redshifts. This explains why the trajectories in the 
three different cases in Figure 2b appear to be nearly identical 
copies of one another, shifted in redshift. 

Note that this prescription only partially captures the clus- 
tering of minihalos. Because of the division into concen- 
tric shells, the radial dimension of the clustering around each 
point is well captured, while the clustering of halos in the 
tangential direction, within each thin shell, is ignored. The 
model, therefore, performs best while ionized bubbles are 
small; otherwise, the effects of clustering are underestimated. 
This means that the predicted minihalo suppression is a lower 
limit, or, equivalently, the predicted bubble masses and ion- 
ized fractions upper limits. 

The probability that a randomly selected point in space will 
land in an ionized bubble of mass M to M+dM is equal to the 
fraction of the total volume (or mass) of the IGM occupied by 
such bubbles, f{M,z) = dx\{M,z), where f(M,z) is the fraction 
of bubble mass histories falling between M and M + dM at 
redshift z- The ionized fraction of the Universe is then 



xi{z)=J2nM,z). 



(11) 



We have performed the above calculation for an ensemble of 
NtofA = 60,000 realizations of the density profile; this implies 
that we can resolve a minimum global ionization fraction of 

^i'-lMotal-lO-^ 

3. RESULTS 

3.1. Universal ionized fraction 

Figure 1 shows the reionization history computed from our 
models. In addition to the unbiased feedback, all-halos, and 
large-halos-only cases discussed previously, the figure also 
displays the result in the case with biased feedback (thick 
solid black curve). Adding clustering dramatically reduces 
the ionized fraction at high redshift. The biased history de- 
parts from the unbiased history before x; = 0.01 and reaches a 
maximum fractional suppression relative to the unbiased case 
of a factor of 11, when i'l.bias = 0.03 vs. Xi.unWas = 0.33 at 



z = 19.2, and a maximum absolute difference of 0.56 when 

^i.bias = 0.22 vs. Xi^unbias = 0.76 at z = 13.6. 

The figure shows that the unbiased and biased reionization 
histories (which have the same efficiency, C = 12) reach per- 
colation at the same redshift; i.e. they both meet at Xi = 1 
at z -^ 11.1. This is an artifact of the modeling procedure, 
and is not physical. It arises because the largest, outermost 
shells in each density profile will always tend to have over- 
densities very close to (5 = 0, so their ionization histories will 
match that of an unbiased Universe, and will therefore be fully 
ionized when the unbiased Universe is fully ionized. In real- 
ity, feedback would always delay percolation, as long as the 
contribution of minihalos to the ionization is not negligible 
(dfcow.mmi/dz ^ d fcoW.Wge / d z; See Fig. 7 below). We also em- 
phasize that our neglect of clustering in the tangential direc- 
tion within each shell will always increase the universal ion- 
ized fraction over what it would be if the clustering in the tan- 
gential direction was also included. Therefore, the difference 
between the unbiased and biased histories is a lower limit, and 
our conclusion that clustering significantly increases minihalo 
suppression is conservative and robust. 

Since observations of Gunn-Peterson troughs in quasar 
spectra constrain the end of reionization to zeor ~ 6, and since 
the value of the efficiency parameter <^ is uncertain, a better 
comparison might be between unbiased and biased trajecto- 
ries with fixed zeor, rather than fixed (^. The artifact men- 
tioned above means that our biased feedback model becomes 
unreliable near percolation, making such an exact compari- 
son impossible. However, we know that our biased feedback 
calculation together with the large-halo-only curve constitute 
upper and lower bounds on the true ionized fraction in the bi- 
ased feedback case. Therefore, if we calculate an unbiased 
feedback history with an efficiency factor lowered (we find, 
to (^ = 5) so that it reaches full ionization at the latest possi- 
ble redshift (zeor ~ 7.7, matching the large-halos-only case), 
then it, together with the (^ = 12 unbiased feedback curve, con- 
stitute commensurable upper and lower bounds. 

In other words, in the context of a fixed zeor, the mini- 
mum possible effect at high redshift of adding clustering is 
shown by the difference between the unbiased feedback curve 
with ( = 5 and the biased feedback curve with ( = 12. The 
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comparison of the thick solid (black) and the light-weight, 
long-dashed (purple) curves in Figure 1 clearly demonstrates 
that minihalo clustering has a significant effect on the high- 
redshift tail of reionization. In Figure lb, which plots the log- 
arithm of the ionized fraction versus redshift, it can be seen 
that the biased feedback curve initially follows the all-halo 
curve, as expected. However, at z ~ 30, the ionized fraction 
reaches a plateau, indicating that nearly all minihalo forma- 
tion is suppressed (see the next subsection); the ionized frac- 
tion does not begin to increase again until the contribution 
from large halos becomes significant. 

3.2. Ionized bubble sizes 
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Fig. 3. — The size distribution of HII regions shown for,fi = 0.01,0.1,0.75 
(left to right in each set). The solid curves include biased feedback, the short- 
dashed curves include all halos and the long-dashed curves are for large halos 
only. The comoving bubble radius is indicated on the bottom axis, and the 
equivalent mass on the top axis. The redshifts corresponding to the ionized 
fractions hsted above are z = 31.7,24.3, 16.3 with all halos, 30.7, 15.3, 11.4 
with biased feedback, and 17.3, 12.9, 8.4 with large halos only. 



In order to better understand the results in the previous sec- 
tion, we here examine the size-distribution of the ionized bub- 
bles. The fraction of points lying in bubbles between mass M 
and M+dM is 



f(M,z) = dx,(z) = V(M)^dM 
dM 



(12) 



where V{M) is the volume of a region of mass M, and 
drtiiii/dM is the number density of ionized bubbles per unit 
mass between M and M+dM. We solve this equation to find 
the mass distribution of the bubbles, given the histogram of 
HII region masses derived from our ensemble of density pro- 
files. Figure 3 plots the size distribution of ionized bubbles, 
normalized by the average ionized fraction of the Universe 
(after FZH), V(M) dnnii/dln(r) x"', for jc, = 0.01,0.1,0.75. 
This function is chosen so that the area under a segment of 
the curve (in log space) is equal to the fraction of ionized gas 
contained in bubbles in that radius (or mass) interval. 

The leftmost set of curves in the figure shows that when 
there is very little ionized gas {xi = 0.01), it is mostly con- 
tained in bubbles around individual halos. Bubble sizes are 
therefore close to Mj = C^min(z)^ the size of a bubble around 
a single, minimum-mass halo. For minihalos. Mi ^ 10'' Mq, 
while for large halos Mi '^ 10** M0. The equivalent comoving 
radii are ^ 10"^ Mpc and ^ 10"' Mpc. 



At the intermediate value ici = 0.1, the feedback case has 
developed a double-peaked distribution. The lower peak cor- 
responds to bubbles close to Mi = (M^iniz) for minihalos, 
while the upper peak corresponds to larger bubbles caused 
by source clustering. Note that this bimodal distribution is 
unique to the biased feedback scenario: since the minihalo 
contribution in ionized bubbles is greatly attenuated, the num- 
ber of isolated bubbles of mass Mi is boosted. By the time 
reionization is nearing completion, xi = 0.75, the lower-radius 
peak has nearly disappeared and the feedback case has a peak 
radius that exceeds the other two cases. It is easy to under- 
stand why the bubble size in the feedback case exceeds that 
of the all-halos case: the contribution of minihalos is atten- 
uated and the larger bubble size reflects the stronger cluster- 
ing of the larger halos. Similar effects emerge when one im- 
poses a mass-dependent ionizing efficiency (Furlanetto et al. 
2006). Note that, in fact, at x, > 0.4, the mean bubble size in 
the feedback case even exceeds that in the large-halos-only 
case. We caution, however, that this result may be unphysical: 
as discussed above, our results become quantitatively unreli- 
able at low redshift. In particular,the mean ionized fraction of 
Xi « 0.4 is reached at the redshift below which the minihalo 
contribution starts to increase (see Fig. 7), suggesting that 
our minihalo suppression becomes only a lower limit below 
these redshifts. It is not clear whether a more accurate treat- 
ment, with more minihalo suppression, would further increase 
or decrease the mean bubble size. Finally, we note that that 
in the biased feedback case, the distribution of bubble sizes 
broadens as reionization proceeds, whereas without feedback 
the bubble size distribution narrows, for fairly general reasons 
(Furlanetto et al. 2006); this may be another signature of the 
unphysical increase in minihalo contribution. 

To gain further intuition, we explore how the mean bubble 
size and the number density of bubbles evolves with redshift 
in each case. We define a characteristic bubble size as the 
mean mass-weighted bubble radius: 



rVir) — ; — dr. 
dr 



(13) 



When the evolution of the characteristic size is considered in 
redshift (Figure 4a), r with feedback is intermediate between 
the all-halo and large-halo-only cases, both in the very be- 
ginning (at z <^ 30, when isolated bubbles dominate), and also 
later (at z ^ 18, when enough large halos have formed to drive 
the distribution). However, when r is considered against the 
ionized fraction xi (Figure 4b), the radius in the feedback case 
quickly exceeds that of both "extremes", although as noted 
above, the mild excess above the large halo only case may not 
be physical. Note that all three curves should approach the 
Hubble radius as xi approaches unity. The unphysical conver- 
gence at M = 10"^Mq occurs because this is the largest mass 
scale included in the numerical modeling. For reference, in 
Figure 4, we also show the correlation lengths, defined by 
b(Mi)b(M2)£,{r) = 1, where b(M) is the Hnear bias of halos 
(Sheth & Tormen 1999), and ^(r) is the correlation function 
of mass. 

The total number density of ionized bubbles, neii = 
J{dnuu/dr)dr, is plotted versus redshift and ionized fraction 
in Figure 5. Comparing this to Figure 4 allows us to describe 
the growth of HII bubbles during the reionization process. 

During the first phase of reionization, most of the ionized 
gas is contained in small bubbles around isolated halos, or 
groups of a few halos. The increase in the bubble size is 
mostly due to the formation of new, larger halos, and the num- 
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Fig. 4. — The mean HII region size versus redshift (a) or versus Universal ionized fraction (b) for the large-halo-only (long-dashed), all-halo (short-dashed), 
and biased feedback (solid) cases. Also shown (nearly-flat curves) is the minimum size A/j for an isolated bubble around a minimum-mass halo for each case. 
The three solid curves that curve downwards at low redshift are (from top to bottom) the correlation lengths for large halos, the cross-correlation length between 
large and minihalos, and the correlation length for small halos. 



ber of HII regions climbs steadily. This phase ends when the 
characteristic bubble size reaches the correlation length of the 
relevant halo population. At this point, the second phase be- 
gins. Bubbles start to merge rapidly, and the growth of the 
characteristic bubble size is now primarily due to those merg- 
ers. The number of HII regions declines steadily as the bubble 
sizes grow larger and larger. Reionization is complete when 
the bubble size becomes formally infinite. 

In the feedback case, the transition to the second phase oc- 
curs at a much smaller ionized fraction than in the other two 
cases (note the early peak in Figure 5b). At z « 24 in the 
all-halo case, the bubble size exceeds the correlation length, 
but this is delayed by minihalo suppression in the feedback 
case until z « 20. By this time there are more bubbles (even 
though the ionized fraction is lower), because more isolated 
halos have had time to form in slightly less-dense regions, 
and these bubbles have not yet merged away. 

As a result, we find that minihalo suppression slows the 
growth of the bubbles and delays their mergers, while allow- 
ing more individual bubbles to form. Because the bubbles are 
smaller, they can be more tightly "packed", i.e. their sepa- 
rations are smaller, so when they do begin to overlap, a very 
rapid period of bubble growth by mergers ensues (z « 22 to 
z ~ 20). Because the mergers occur so rapidly (and because 
of high minihalo suppression), the ionized fraction grows very 
little, while the bubble size grows by over an order of magni- 
tude (note the sharpness of the peak in Figure 5b and the rapid 
rise in radius in Figure 4b). This means that the mean bubble 
size for the feedback case ends up very high, for any given 
ionized fraction at x; ^ 0.05. 

3.3. Mean bias factor 

Our main result above can be summarized as follows: a 
much larger fraction of minihalos are suppressed when their 
clustering is included, relative to the value of Xi expected if the 
minihalos were randomly distributed. Here we compare the 
"boost" in the fraction of suppressed halos we find to naive 
expectations, using a linear halo-bias model. 

Let us first assume that the average ionized fraction of the 
Universe evolves according to this differential equation. 



dxi 
dz 



dx\ 



i, large 



+ {i-B{z)x-; 



dxii 



(14) 



dz dz 

where B{z) is the effective "boost", i.e. the enhancement in 



the number of minihalos forming in ionized regions, relative 
to Xi, due to their biased distribution. Since we know each 
dxi/dz and x; from our simulation runs, this equation serves 
as a definition of B{z). Figure 6 plots this factor, calculated 
from the biased feedback ionization history using 



5(z)=^(l- 

Xi 



dxi/dz-dxu^igf./dz 



(15) 



^■^i.mim/ aZ 

Because the dx-Jdz > dxi^xsign/dz, we know that B > x~' . Note 
that B = xj"' represents the limiting case in which all minihalos 
are suppressed. 

Figure 6 shows that at the very beginning of reionization, 
the enhancement is relatively low (much less than iij"') since 
few bubbles have grown large enough to encompass more 
than one halo. However, the boost factor quickly climbs to 



Bfe) 



50, and at z ^ 28, nearly all minihalos are sup- 



pressed. This corresponds to the plateau of the solid curve 
in Figure lb, and the epoch of rapid mergers discussed in 
the previous subsection. The enhancement subsequently de- 
creases as the ionized fraction increases and the ionized re- 
gions are sampling regions whose density is closer to the 
global average. Values of B{z) < 1 are artifacts of the model, 
caused by the unphy sically rapid end of reionization discussed 
in §3.1 above. 

Naively, we might expect that the enhancement factor 
would be roughly equal to the enhancement in the average 
number of halos within a bubble radius of an average halo. 



B 



neighbor ~ 



?>b{Mx)b{M2) 



r Cn(r)dr, 



(16) 



where b(M) is the bias of halos of mass M, ^,„(r) is the mass 
correlation function, and r^in is a minimum radius chosen to 
exclude overlapping halos (though the calculation is quite in- 
sensitive to its exact value). Figure 6 shows this enhance- 
ment for large halos near large halos, large halos near miniha- 
los, and minihalos near minihalos (using the minimum masses 
of each halo type). At the earliest stages of reionization the 
lowest curve — the enhancement in the number of minihalos 
within a bubble radius of an average minihalo — should be 
directly applicable. The figure, however, shows that this un- 
derestimates the actual bias enhancement by nearly an order 
of magnitude. 

Several effects may contribute to this interesting disparity. 
First, we are interested in small scales (the ionized mass is 
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Fig. 5. — The comoving number density of HII regions versus redshift (a) or versus Universal ionized fraction (b) for the large-halo-only (long-dashed), 
all-halo (short-dashed), and feedback (solid) cases. The plots are normalized by the maximum number density for each case, which is 1.0 X Iff* Mpc"^ with all 
halos, 1.5 X lO** Mpc"^ with biased feedback, and 2.8 Mpc"' with large halos only. 



only 12 times the halo mass), where using the linear power 
spectrum, as we did in equation (16), will underpredict the 
clustering. As shown by Iliev et al. (2003), the effects of 
non-linear bias can boost the clustering of minihalos at high 
redshift on small scales substantially over the standard lin- 
ear result (for instance, at z ^ 20, when our typical bubble 
size is < O.lMpc, Fig. 2 of Iliev et al. (2003) shows that 
non-linear bias can exceed the standard linear result by an 
order of magnitude). They show that a non-linear analytic 
technique based on extended Press-Schechter theory provides 
an accurate match to numerical simulations; our excursion- 
set based formalism, although different in detail, should like- 
wise capture non-linear bias. Such effects could in princi- 
ple be largely responsible for the difference between our re- 
sults and equation (16), although the overall contribution is 
difficult to assess since volume-exclusion effects will reduce 
clustering on small scales and will compensate somewhat for 
our neglecting nonlinear effects. It clearly dominates at the 
beginning of reionization, when most bubbles have a typical 
size M] = C-'^min(z)- Sccond, equation (16) gives the increase 
in the number of neighbors centered around a single halo, 
whereas in our case, the ionized region is already known to 
contain several halos. Third, equation (16) averages the bias 
over all halos at a given redshift; in our case, we are interested 
in the number of new halos forming around a set of halos that 
had formed earlier. Fourth, the bubble-size distribution in the 
feedback case always has a tail extending to small sizes, so 
using the mass weighted mean radius r will underestimate the 
bias; in principle, one should integrate over the distribution of 
bubble sizes. 

In summary, we conclude that given all the caveats 
and complications mentioned above, either a tailored semi- 
analytic approach or a naive estimate based on linear the- 
ory will both give misleading underestimates of the impact 
of clustering. 

3.4. Minihalo suppression and electron scattering optical 

depth 

In Figure 7, we show the fraction of the reionization con- 
tributed by minihalos as a function of redshift in the biased 
and unbiased feedback cases. For the reasons discussed pre- 
viously, the minihalo fraction calculated here for the biased 
feedback case is an upper Umit. In particular, the apparent up- 



TABLE 1 
Optical depth and end-of-reionization redshift for 

EACH scenario 



Scenario: 


all 


unbiased 


unbiased 


clustered 


large only 


c 


12 


12 


5 


12 


12 


r 


0.21 


0.20 


0.14 


0.12 


0.08 


ZEOR 


15.1 


11.1 


7.8 


11.1 


7.7 



turn of the minihalo fraction at low redshift is unphysical. The 
most conservative conclusion is that less than 60% of the full 
ionization would be contributed by minihalos in a more exact 
treatment of clustered feedback. It is likely, however, that in 
the absence of the unphysical upturn, the minihalo contribu- 
tion would fall significantly below this value at z <^ 15. 

Finally, in light of the recent measurement in the three-year 
WMAP data, it is interesting to compute the impact of mini- 
halo suppression on the optical depth t to electron scattering. 
We compute the value of r using the evolution of the ion- 
ized fraction of hydrogen in our models. We ignore additional 
electrons from ionized Helium (if the fraction of singly ion- 
ized helium would track that of hydrogen, r would increase 
by - 8%). 

The optical depths for each history are reported in Table 1 . 
We naturally find the earliest ionization in the no feedback 
(all-halos) case, so this limiting model has the highest opti- 
cal depth. With unbiased feedback, reionization is delayed, 
and we find a slightly lower optical depth. The lower effi- 
ciency unbiased case has even lower optical depth. And our 
history with clustered feedback results in lower optical depth 
than any other scenario, except the case in which we elimi- 
nate all minihalos. Note that since the biased feedback case 
is an upper limit on the ionized fraction and reaches Xi = 1 
artificially early, the calculated optical depth is also an upper 
limit in this case. On the other hand, the optical depth in the 
large-halo-only case constitutes a lower limit to the true bi- 
ased feedback case. 

In conclusion. Table 1 shows that clustering reduces the op- 
tical depth from somewhere in the range of 0.14 < r < 0.20 
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to somewhere in the range of 0.08 < t < 0.12. 

4. DISCUSSION 
4.1. Caveats 

Obviously many simplifications must be made to construct 
a semi-analytical model of a process as complicated as reion- 
ization. We have already discussed the fact that our results 
for feedback with clustering constitute an upper limit on the 
ionized fraction due to the averaging of the suppression effect 
over spherical shells and the dominance of large shells when 
the bubbles are large. However, we only start to see obviously 
unphysical effects (B < 1, or an increasing minihalo contribu- 
tion to reionization) at xi > 0.4 or z < 12. We expect that a 
more exact treatment that incorporates feedback with a three- 
dimensional treatment of clustering would stay close to this 
upper limit prior to this redshift. 

Another important oversimplification is that we have ig- 
nored recombinations (except in a very rough way in the ( 
factor). In particular, we ignore the scale dependence of the 
recombination rate; as the bubbles get bigger and the mean 
free path of ionizing photons increases, the gas in the IGM is 
likely ionized up to larger overdensities, and the importance 
of recombinations will increase. The omission of recombina- 
tion may also be important in the "stalling" regime we find 
occurs at high redshift z ^ 20: since there are few ionizing 
photons being produced and the ionized bubbles are in rela- 
tively dense regions, recombinations would likely lower the 
ionized fraction. Furlanetto & Oh (2005) describe a method 
for including recombinations in the FZH formalism, which 
might be adapted for use in this model. Although their results 
suggest that recombinations are only important in the very late 
stages of reionization, the larger bubble sizes produced in the 
feedback case may make recombinations important somewhat 
earlier. However, since the late-time behavior of our biased 
feedback model when x, ^ 1 is in any case unreliable (as dis- 
cussed in §2), we have foregone this extra complication in the 
present paper. 

Furlanetto et al. (2006) and Cohn & Chang (2006) have ex- 
plored several generalizations of the FZH formalism, includ- 
ing the use of a more suitable halo mass function, stochastic 
fluctuations in the halo distribution, a mass-dependent effi- 
ciency factor C(M), and the history of halo mergers within 
ionized bubbles. Their results suggest that the first two mod- 
ifications would have relatively minor effects on our results. 
A mass-dependent efficiency, or associating sources with halo 
mergers, however, could both change the importance of mini- 
halo suppression, since it could change the relative contribu- 
tions of minihalos and large halos to reionization. In particu- 
lar, if the efficiency of ionizing photon production in miniha- 
los were reduced enough (by internal feedback, for instance), 
minihalo suppression in ionized regions would become unim- 
portant. 

Our specific assumptions about minihalo suppression are 
also somewhat idealized. For instance, there will clearly not 
be a sharp cut-off of minihalo formation at precisely 10^ K. 



There is also some uncertainty about the exact value of the 
minimum virial temperature, which may have some environ- 
mental dependence, in any case. Changing the dividing line 
between minihalos and large halos will change their relative 
populations and affect the precise reionization history, but 
small changes should not affect the importance of minihalo 
suppression as a feedback mechanism in general. 

4.2. Possible extensions 

There are several extensions to the formalism presented 
here that might be profitably explored in the future, such as 
including a mass-dependent efficiency factor, C(M), and in 
particular using different efficiencies for large and minihalos. 
There are possibilities for improving the performance of the 
model at high xi, which would allow us to better constrain the 
end-of-reionization redshift and the optical depth, or, more 
importantly, to calculate r under various assumptions given a 
fixed zeor- As mentioned above, our formalism can be eas- 
ily adapted to treat a variety of feedback mechanisms other 
than minihalo suppression; our results suggest that the im- 
pact of these other feedback mechanisms will generically be 
increased by the presence of clustering, as well. We plan to 
apply the formalism discussed here to address the enrichment 
of the IGM with metals in a future paper. 

5. CONCLUSIONS 

Using a unique semi-analytical technique, we have demon- 
strated that minihalo suppression in ionized regions is en- 
hanced by clustering, and that such enhanced suppression 
substantially reduces the high-redshift tail of reionization. 
Other feedback mechanisms may be similarly enhanced by 
clustering, since feedback process tend to operates over a lim- 
ited distance, especially when sources are young. 

The compressed ionization history produced by minihalo 
suppression reduces the optical depth r compared to cases 
without clustering or without feedback. The reduction in r 
occurs both with fixed ionization efficiency factor C, and with 
fixed end-of-reionization redshift zeor- 

Haiman & Bryan (2006) have shown that substantial mini- 
halo suppression is necessary to simultaneously meet the con- 
straints imposed by Gunn-Peterson troughs on zeor and by 
the value of r measured in the three-year data of the WMAP 
experiment. We conclude that a simple feedback mechanism 
— such as preventing the formation of minihalos in ionized 
regions — can produce such suppression when the highly- 
clustered distribution of these sources is taken into account. 
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